Stokes I Simultaneous Image and Instrument Modeling

In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 by simultaneously creating an image and model for the instrument. By instrument model, we mean something akin to self-calibration in traditional VLBI imaging terminology. However, unlike traditional self-cal, we will at each point in our parameter space effectively explore the possible self-cal solutions. This will allow us to constrain and marginalize over the instrument effects, such as time variable gains.

To get started we load Comrade.

using Comrade
  Activating project at `~/work/Comrade.jl/Comrade.jl/examples`
using Pyehtim
using LinearAlgebra

For reproducibility we use a stable random number genreator

using StableRNGs
rng = StableRNG(123)
StableRNGs.LehmerRNG(state=0x000000000000000000000000000000f7)

Load the Data

To download the data visit https://doi.org/10.25739/g85n-f134 First we will load our data:

obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "examples", "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
Python: <ehtim.obsdata.Obsdata object at 0x7fcb6a7c84c0>

Now we do some minor preprocessing:

  • Scan average the data since the data have been preprocessed so that the gain phases coherent.
  • Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.
obs = scan_average(obs.add_fractional_noise(0.02))
Python: <ehtim.obsdata.Obsdata object at 0x7fcb89242ec0>

Now we extract our complex visibilities.

dvis = extract_table(obs, ComplexVisibilities())
EHTObservation{Float64,Comrade.EHTVisibilityDatum{Float64}, ...}
  source: M87
  mjd: 57849
  frequency: 2.29070703125e11
  bandwidth: 1.856e9
  stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples: 284

##Building the Model/Posterior

Now, we must build our intensity/visibility model. That is, the model that takes in a named tuple of parameters and perhaps some metadata required to construct the model. For our model, we will use a raster or ContinuousImage for our image model. The model is given below:

function sky(θ, metadata)
    (;fg, c, σimg) = θ
    (;ftot, K, meanpr, grid, cache) = metadata
    # Transform to the log-ratio pixel fluxes
    cp = meanpr .+ σimg.*c.params
    # Transform to image space
    rast = (ftot*(1-fg))*K(to_simplex(CenteredLR(), cp))
    img = IntensityMap(rast, grid)
    m = ContinuousImage(img, cache)
    # Add a large-scale gaussian to deal with the over-resolved mas flux
    g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(ftot*fg))
    return m + g
end
sky (generic function with 1 method)

Unlike other imaging examples (e.g., Imaging a Black Hole using only Closure Quantities) we also need to include a model for the instrument, i.e., gains as well. The gains will be broken into two components

  • Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
  • Gain phases which are more difficult to constrain and can shift rapidly.
function instrument(θ, metadata)
    (; lgamp, gphase) = θ
    (; gcache, gcachep) = metadata
    # Now form our instrument model
    gvis = exp.(lgamp)
    gphase = exp.(1im.*gphase)
    jgamp = jonesStokes(gvis, gcache)
    jgphase = jonesStokes(gphase, gcachep)
    return JonesModel(jgamp*jgphase)
end
instrument (generic function with 1 method)

The model construction is very similar to Imaging a Black Hole using only Closure Quantities, except we include a large scale gaussian since we want to model the zero baselines. For more information about the image model please read the closure-only example. Let's discuss the instrument model Comrade.JonesModel. Thanks to the EHT pre-calibration, the gains are stable over scans. Therefore, we can model the gains on a scan-by-scan basis. To form the instrument model, we need our

  1. Our (log) gain amplitudes and phases are given below by lgamp and gphase
  2. Our function or cache that maps the gains from a list to the stations they impact gcache.
  3. The set of Comrade.JonesPairs produced by jonesStokes

These three ingredients then specify our instrument model. The instrument model can then be combined with our image model cimg to form the total JonesModel.

Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of view. Typically 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.

npix = 32
fovx = μas2rad(150.0)
fovy = μas2rad(150.0)
7.27220521664304e-10

Now let's form our cache's. First, we have our usual image cache which is needed to numerically compute the visibilities.

grid = imagepixels(fovx, fovy, npix, npix)
buffer = IntensityMap(zeros(npix, npix), grid)
cache = create_cache(NFFTAlg(dvis), buffer, BSplinePulse{3}())
VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im  …  0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Second, we now construct our instrument model cache. This tells us how to map from the gains to the model visibilities. However, to construct this map, we also need to specify the observation segmentation over which we expect the gains to change. This is specified in the second argument to jonescache, and currently, there are two options

  • FixedSeg(val): Fixes the corruption to the value val for all time. This is usefule for reference stations
  • ScanSeg(): which forces the corruptions to only change from scan-to-scan
  • TrackSeg(): which forces the corruptions to be constant over a night's observation

For this work, we use the scan segmentation for the gain amplitudes since that is roughly the timescale we expect them to vary. For the phases we use a station specific scheme where we set AA to be fixed to unit gain because it will function as a reference station.

gcache = jonescache(dvis, ScanSeg())
gcachep = jonescache(dvis, ScanSeg(); autoref=SEFDReference((complex(1.0))))

using VLBIImagePriors

Now we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior Since we are using a Gaussian Markov random field prior we need to first specify our mean image. This behaves somewhat similary to a entropy regularizer in that it will start with an initial guess for the image structure. For this tutorial we will use a a symmetric Gaussian with a FWHM of 60 μas

fwhmfac = 2*sqrt(2*log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0)./fwhmfac))
imgpr = intensitymap(mpr, grid)
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 NamedDimsArray(::Matrix{Float64}, (:X, :Y)):
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8
 (-3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
 (-3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
 (-2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
 (-2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
 (-2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46586e-6        5.12225e-6     2.46586e-6
  (2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
  (2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
  (2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
  (3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
  (3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
  (3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8

Now since we are actually modeling our image on the simplex we need to ensure that our mean image has unit flux

imgpr ./= flux(imgpr)
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 Matrix{Float64}:
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.38049e-8        1.3254e-7      6.38049e-8
 (-3.29522e-10)     1.3254e-7         2.75321e-7     1.3254e-7
 (-3.06796e-10)     2.62224e-7        5.4471e-7      2.62224e-7
 (-2.84071e-10)     4.94121e-7        1.02642e-6     4.94121e-7
 (-2.61345e-10)     8.86803e-7   …    1.84213e-6     8.86803e-7
 (-2.38619e-10)     1.51585e-6        3.14882e-6     1.51585e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46784e-6        5.12636e-6     2.46784e-6
  (2.38619e-10)     1.51585e-6        3.14882e-6     1.51585e-6
  (2.61345e-10)     8.86803e-7   …    1.84213e-6     8.86803e-7
  (2.84071e-10)     4.94121e-7        1.02642e-6     4.94121e-7
  (3.06796e-10)     2.62224e-7        5.4471e-7      2.62224e-7
  (3.29522e-10)     1.3254e-7         2.75321e-7     1.3254e-7
  (3.52247e-10)     6.38049e-8        1.3254e-7      6.38049e-8

and since our prior is not on the simplex we need to convert it to unconstrained or real space.

meanpr = to_real(CenteredLR(), Comrade.baseimage(imgpr))
32×32 Matrix{Float64}:
 -7.55422  -6.82317  -6.14085  -5.50727   …  -6.14085  -6.82317  -7.55422
 -6.82317  -6.09211  -5.4098   -4.77622      -5.4098   -6.09211  -6.82317
 -6.14085  -5.4098   -4.72748  -4.0939       -4.72748  -5.4098   -6.14085
 -5.50727  -4.77622  -4.0939   -3.46032      -4.0939   -4.77622  -5.50727
 -4.92243  -4.19137  -3.50906  -2.87548      -3.50906  -4.19137  -4.92243
 -4.38632  -3.65527  -2.97295  -2.33937   …  -2.97295  -3.65527  -4.38632
 -3.89895  -3.1679   -2.48558  -1.852        -2.48558  -3.1679   -3.89895
 -3.46032  -2.72927  -2.04695  -1.41337      -2.04695  -2.72927  -3.46032
 -3.07043  -2.33937  -1.65705  -1.02348      -1.65705  -2.33937  -3.07043
 -2.72927  -1.99821  -1.3159   -0.682317     -1.3159   -1.99821  -2.72927
  ⋮                                       ⋱             ⋮        
 -3.07043  -2.33937  -1.65705  -1.02348      -1.65705  -2.33937  -3.07043
 -3.46032  -2.72927  -2.04695  -1.41337      -2.04695  -2.72927  -3.46032
 -3.89895  -3.1679   -2.48558  -1.852     …  -2.48558  -3.1679   -3.89895
 -4.38632  -3.65527  -2.97295  -2.33937      -2.97295  -3.65527  -4.38632
 -4.92243  -4.19137  -3.50906  -2.87548      -3.50906  -4.19137  -4.92243
 -5.50727  -4.77622  -4.0939   -3.46032      -4.0939   -4.77622  -5.50727
 -6.14085  -5.4098   -4.72748  -4.0939       -4.72748  -5.4098   -6.14085
 -6.82317  -6.09211  -5.4098   -4.77622   …  -5.4098   -6.09211  -6.82317
 -7.55422  -6.82317  -6.14085  -5.50727      -6.14085  -6.82317  -7.55422

Now we can form our metadata we need to fully define our model.

metadata = (;ftot=1.1, K=CenterImage(imgpr), meanpr, grid, cache, gcache, gcachep)
(ftot = 1.1, K = VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}([0.9944957386363662 -0.005326704545454548 … 0.005326704545454572 0.005504261363636381; -0.005326704545454548 0.9948393969941349 … 0.005160603005865089 0.005326704545454565; … ; 0.005326704545454572 0.005160603005865089 … 0.9948393969941348 -0.005326704545454546; 0.005504261363636381 0.005326704545454565 … -0.005326704545454546 0.9944957386363638], (32, 32)), meanpr = [-7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; … ; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; -7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378], grid = GriddedKeys{(:X, :Y)}
	X: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
	Y: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
, cache = VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im  …  0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]), gcache = JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25)), gcachep = JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))]))

We will also fix the total flux to be the observed value 1.1. This is because total flux is degenerate with a global shift in the gain amplitudes making the problem degenerate. To fix this we use the observed total flux as our value.

Moving onto our prior, we first focus on the instrument model priors. Each station requires its own prior on both the amplitudes and phases. For the amplitudes we assume that the gains are apriori well calibrated around unit gains (or 0 log gain amplitudes) which corresponds to no instrument corruption. The gain dispersion is then set to 10% for all stations except LMT, representing that we expect 10% deviations from scan-to-scan. For LMT we let the prior expand to 100% due to the known pointing issues LMT had in 2017.

using Distributions
using DistributionsAD
distamp = station_tuple(dvis, Normal(0.0, 0.1); LM = Normal(1.0))
(AA = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AP = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AZ = Distributions.Normal{Float64}(μ=0.0, σ=0.1), JC = Distributions.Normal{Float64}(μ=0.0, σ=0.1), LM = Distributions.Normal{Float64}(μ=1.0, σ=1.0), PV = Distributions.Normal{Float64}(μ=0.0, σ=0.1), SM = Distributions.Normal{Float64}(μ=0.0, σ=0.1))

For the phases, as mentioned above, we will use a segmented gain prior. This means that rather than the parameters being directly the gains, we fit the first gain for each site, and then the other parameters are the segmented gains compared to the previous time. To model this we break the gain phase prior into two parts. The first is the prior for the first observing timestamp of each site, distphase0, and the second is the prior for segmented gain ϵₜ from time i to i+1, given by distphase. For the EHT, we are dealing with pre-2*rand(rng, ndim) .- 1.5calibrated data, so often, the gain phase jumps from scan to scan are minor. As such, we can put a more informative prior on distphase.

Warning

We use AA (ALMA) as a reference station so we do not have to specify a gain prior for it.

distphase = station_tuple(dvis, DiagonalVonMises(0.0, inv(π^2)))
(AA = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AP = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AZ = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), JC = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), LM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), PV = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), SM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688))

In addition we want a reasonable guess for what the resolution of our image should be. For radio astronomy this is given by roughly the longest baseline in the image. To put this into pixel space we then divide by the pixel size.

beam = beamsize(dvis)
rat = (beam/(step(grid.X)))
5.279832635689346

To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. This drastically improves the usual N^3 scaling you get from usual Gaussian Processes.

crcache = MarkovRandomFieldCache(meanpr)
VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}(sparse([1, 2, 32, 33, 993, 1, 2, 3, 34, 994  …  31, 991, 1022, 1023, 1024, 32, 992, 993, 1023, 1024], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, 119, 119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, 121, 122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, 124, 124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, 126, 126, 127, 127, 127, 127, 127, 128, 128, 128, 128, 128, 129, 129, 129, 129, 129, 130, 130, 130, 130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 132, 132, 133, 133, 133, 133, 133, 134, 134, 134, 134, 134, 135, 135, 135, 135, 135, 136, 136, 136, 136, 136, 137, 137, 137, 137, 137, 138, 138, 138, 138, 138, 139, 139, 139, 139, 139, 140, 140, 140, 140, 140, 141, 141, 141, 141, 141, 142, 142, 142, 142, 142, 143, 143, 143, 143, 143, 144, 144, 144, 144, 144, 145, 145, 145, 145, 145, 146, 146, 146, 146, 146, 147, 147, 147, 147, 147, 148, 148, 148, 148, 148, 149, 149, 149, 149, 149, 150, 150, 150, 150, 150, 151, 151, 151, 151, 151, 152, 152, 152, 152, 152, 153, 153, 153, 153, 153, 154, 154, 154, 154, 154, 155, 155, 155, 155, 155, 156, 156, 156, 156, 156, 157, 157, 157, 157, 157, 158, 158, 158, 158, 158, 159, 159, 159, 159, 159, 160, 160, 160, 160, 160, 161, 161, 161, 161, 161, 162, 162, 162, 162, 162, 163, 163, 163, 163, 163, 164, 164, 164, 164, 164, 165, 165, 165, 165, 165, 166, 166, 166, 166, 166, 167, 167, 167, 167, 167, 168, 168, 168, 168, 168, 169, 169, 169, 169, 169, 170, 170, 170, 170, 170, 171, 171, 171, 171, 171, 172, 172, 172, 172, 172, 173, 173, 173, 173, 173, 174, 174, 174, 174, 174, 175, 175, 175, 175, 175, 176, 176, 176, 176, 176, 177, 177, 177, 177, 177, 178, 178, 178, 178, 178, 179, 179, 179, 179, 179, 180, 180, 180, 180, 180, 181, 181, 181, 181, 181, 182, 182, 182, 182, 182, 183, 183, 183, 183, 183, 184, 184, 184, 184, 184, 185, 185, 185, 185, 185, 186, 186, 186, 186, 186, 187, 187, 187, 187, 187, 188, 188, 188, 188, 188, 189, 189, 189, 189, 189, 190, 190, 190, 190, 190, 191, 191, 191, 191, 191, 192, 192, 192, 192, 192, 193, 193, 193, 193, 193, 194, 194, 194, 194, 194, 195, 195, 195, 195, 195, 196, 196, 196, 196, 196, 197, 197, 197, 197, 197, 198, 198, 198, 198, 198, 199, 199, 199, 199, 199, 200, 200, 200, 200, 200, 201, 201, 201, 201, 201, 202, 202, 202, 202, 202, 203, 203, 203, 203, 203, 204, 204, 204, 204, 204, 205, 205, 205, 205, 205, 206, 206, 206, 206, 206, 207, 207, 207, 207, 207, 208, 208, 208, 208, 208, 209, 209, 209, 209, 209, 210, 210, 210, 210, 210, 211, 211, 211, 211, 211, 212, 212, 212, 212, 212, 213, 213, 213, 213, 213, 214, 214, 214, 214, 214, 215, 215, 215, 215, 215, 216, 216, 216, 216, 216, 217, 217, 217, 217, 217, 218, 218, 218, 218, 218, 219, 219, 219, 219, 219, 220, 220, 220, 220, 220, 221, 221, 221, 221, 221, 222, 222, 222, 222, 222, 223, 223, 223, 223, 223, 224, 224, 224, 224, 224, 225, 225, 225, 225, 225, 226, 226, 226, 226, 226, 227, 227, 227, 227, 227, 228, 228, 228, 228, 228, 229, 229, 229, 229, 229, 230, 230, 230, 230, 230, 231, 231, 231, 231, 231, 232, 232, 232, 232, 232, 233, 233, 233, 233, 233, 234, 234, 234, 234, 234, 235, 235, 235, 235, 235, 236, 236, 236, 236, 236, 237, 237, 237, 237, 237, 238, 238, 238, 238, 238, 239, 239, 239, 239, 239, 240, 240, 240, 240, 240, 241, 241, 241, 241, 241, 242, 242, 242, 242, 242, 243, 243, 243, 243, 243, 244, 244, 244, 244, 244, 245, 245, 245, 245, 245, 246, 246, 246, 246, 246, 247, 247, 247, 247, 247, 248, 248, 248, 248, 248, 249, 249, 249, 249, 249, 250, 250, 250, 250, 250, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 253, 253, 253, 253, 253, 254, 254, 254, 254, 254, 255, 255, 255, 255, 255, 256, 256, 256, 256, 256, 257, 257, 257, 257, 257, 258, 258, 258, 258, 258, 259, 259, 259, 259, 259, 260, 260, 260, 260, 260, 261, 261, 261, 261, 261, 262, 262, 262, 262, 262, 263, 263, 263, 263, 263, 264, 264, 264, 264, 264, 265, 265, 265, 265, 265, 266, 266, 266, 266, 266, 267, 267, 267, 267, 267, 268, 268, 268, 268, 268, 269, 269, 269, 269, 269, 270, 270, 270, 270, 270, 271, 271, 271, 271, 271, 272, 272, 272, 272, 272, 273, 273, 273, 273, 273, 274, 274, 274, 274, 274, 275, 275, 275, 275, 275, 276, 276, 276, 276, 276, 277, 277, 277, 277, 277, 278, 278, 278, 278, 278, 279, 279, 279, 279, 279, 280, 280, 280, 280, 280, 281, 281, 281, 281, 281, 282, 282, 282, 282, 282, 283, 283, 283, 283, 283, 284, 284, 284, 284, 284, 285, 285, 285, 285, 285, 286, 286, 286, 286, 286, 287, 287, 287, 287, 287, 288, 288, 288, 288, 288, 289, 289, 289, 289, 289, 290, 290, 290, 290, 290, 291, 291, 291, 291, 291, 292, 292, 292, 292, 292, 293, 293, 293, 293, 293, 294, 294, 294, 294, 294, 295, 295, 295, 295, 295, 296, 296, 296, 296, 296, 297, 297, 297, 297, 297, 298, 298, 298, 298, 298, 299, 299, 299, 299, 299, 300, 300, 300, 300, 300, 301, 301, 301, 301, 301, 302, 302, 302, 302, 302, 303, 303, 303, 303, 303, 304, 304, 304, 304, 304, 305, 305, 305, 305, 305, 306, 306, 306, 306, 306, 307, 307, 307, 307, 307, 308, 308, 308, 308, 308, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311, 311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313, 313, 314, 314, 314, 314, 314, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 317, 317, 317, 317, 317, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 320, 320, 320, 320, 320, 321, 321, 321, 321, 321, 322, 322, 322, 322, 322, 323, 323, 323, 323, 323, 324, 324, 324, 324, 324, 325, 325, 325, 325, 325, 326, 326, 326, 326, 326, 327, 327, 327, 327, 327, 328, 328, 328, 328, 328, 329, 329, 329, 329, 329, 330, 330, 330, 330, 330, 331, 331, 331, 331, 331, 332, 332, 332, 332, 332, 333, 333, 333, 333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 336, 336, 336, 336, 336, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 339, 339, 339, 339, 339, 340, 340, 340, 340, 340, 341, 341, 341, 341, 341, 342, 342, 342, 342, 342, 343, 343, 343, 343, 343, 344, 344, 344, 344, 344, 345, 345, 345, 345, 345, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 352, 352, 352, 352, 352, 353, 353, 353, 353, 353, 354, 354, 354, 354, 354, 355, 355, 355, 355, 355, 356, 356, 356, 356, 356, 357, 357, 357, 357, 357, 358, 358, 358, 358, 358, 359, 359, 359, 359, 359, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 362, 362, 362, 362, 362, 363, 363, 363, 363, 363, 364, 364, 364, 364, 364, 365, 365, 365, 365, 365, 366, 366, 366, 366, 366, 367, 367, 367, 367, 367, 368, 368, 368, 368, 368, 369, 369, 369, 369, 369, 370, 370, 370, 370, 370, 371, 371, 371, 371, 371, 372, 372, 372, 372, 372, 373, 373, 373, 373, 373, 374, 374, 374, 374, 374, 375, 375, 375, 375, 375, 376, 376, 376, 376, 376, 377, 377, 377, 377, 377, 378, 378, 378, 378, 378, 379, 379, 379, 379, 379, 380, 380, 380, 380, 380, 381, 381, 381, 381, 381, 382, 382, 382, 382, 382, 383, 383, 383, 383, 383, 384, 384, 384, 384, 384, 385, 385, 385, 385, 385, 386, 386, 386, 386, 386, 387, 387, 387, 387, 387, 388, 388, 388, 388, 388, 389, 389, 389, 389, 389, 390, 390, 390, 390, 390, 391, 391, 391, 391, 391, 392, 392, 392, 392, 392, 393, 393, 393, 393, 393, 394, 394, 394, 394, 394, 395, 395, 395, 395, 395, 396, 396, 396, 396, 396, 397, 397, 397, 397, 397, 398, 398, 398, 398, 398, 399, 399, 399, 399, 399, 400, 400, 400, 400, 400, 401, 401, 401, 401, 401, 402, 402, 402, 402, 402, 403, 403, 403, 403, 403, 404, 404, 404, 404, 404, 405, 405, 405, 405, 405, 406, 406, 406, 406, 406, 407, 407, 407, 407, 407, 408, 408, 408, 408, 408, 409, 409, 409, 409, 409, 410, 410, 410, 410, 410, 411, 411, 411, 411, 411, 412, 412, 412, 412, 412, 413, 413, 413, 413, 413, 414, 414, 414, 414, 414, 415, 415, 415, 415, 415, 416, 416, 416, 416, 416, 417, 417, 417, 417, 417, 418, 418, 418, 418, 418, 419, 419, 419, 419, 419, 420, 420, 420, 420, 420, 421, 421, 421, 421, 421, 422, 422, 422, 422, 422, 423, 423, 423, 423, 423, 424, 424, 424, 424, 424, 425, 425, 425, 425, 425, 426, 426, 426, 426, 426, 427, 427, 427, 427, 427, 428, 428, 428, 428, 428, 429, 429, 429, 429, 429, 430, 430, 430, 430, 430, 431, 431, 431, 431, 431, 432, 432, 432, 432, 432, 433, 433, 433, 433, 433, 434, 434, 434, 434, 434, 435, 435, 435, 435, 435, 436, 436, 436, 436, 436, 437, 437, 437, 437, 437, 438, 438, 438, 438, 438, 439, 439, 439, 439, 439, 440, 440, 440, 440, 440, 441, 441, 441, 441, 441, 442, 442, 442, 442, 442, 443, 443, 443, 443, 443, 444, 444, 444, 444, 444, 445, 445, 445, 445, 445, 446, 446, 446, 446, 446, 447, 447, 447, 447, 447, 448, 448, 448, 448, 448, 449, 449, 449, 449, 449, 450, 450, 450, 450, 450, 451, 451, 451, 451, 451, 452, 452, 452, 452, 452, 453, 453, 453, 453, 453, 454, 454, 454, 454, 454, 455, 455, 455, 455, 455, 456, 456, 456, 456, 456, 457, 457, 457, 457, 457, 458, 458, 458, 458, 458, 459, 459, 459, 459, 459, 460, 460, 460, 460, 460, 461, 461, 461, 461, 461, 462, 462, 462, 462, 462, 463, 463, 463, 463, 463, 464, 464, 464, 464, 464, 465, 465, 465, 465, 465, 466, 466, 466, 466, 466, 467, 467, 467, 467, 467, 468, 468, 468, 468, 468, 469, 469, 469, 469, 469, 470, 470, 470, 470, 470, 471, 471, 471, 471, 471, 472, 472, 472, 472, 472, 473, 473, 473, 473, 473, 474, 474, 474, 474, 474, 475, 475, 475, 475, 475, 476, 476, 476, 476, 476, 477, 477, 477, 477, 477, 478, 478, 478, 478, 478, 479, 479, 479, 479, 479, 480, 480, 480, 480, 480, 481, 481, 481, 481, 481, 482, 482, 482, 482, 482, 483, 483, 483, 483, 483, 484, 484, 484, 484, 484, 485, 485, 485, 485, 485, 486, 486, 486, 486, 486, 487, 487, 487, 487, 487, 488, 488, 488, 488, 488, 489, 489, 489, 489, 489, 490, 490, 490, 490, 490, 491, 491, 491, 491, 491, 492, 492, 492, 492, 492, 493, 493, 493, 493, 493, 494, 494, 494, 494, 494, 495, 495, 495, 495, 495, 496, 496, 496, 496, 496, 497, 497, 497, 497, 497, 498, 498, 498, 498, 498, 499, 499, 499, 499, 499, 500, 500, 500, 500, 500, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 503, 503, 503, 503, 503, 504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 506, 506, 506, 506, 506, 507, 507, 507, 507, 507, 508, 508, 508, 508, 508, 509, 509, 509, 509, 509, 510, 510, 510, 510, 510, 511, 511, 511, 511, 511, 512, 512, 512, 512, 512, 513, 513, 513, 513, 513, 514, 514, 514, 514, 514, 515, 515, 515, 515, 515, 516, 516, 516, 516, 516, 517, 517, 517, 517, 517, 518, 518, 518, 518, 518, 519, 519, 519, 519, 519, 520, 520, 520, 520, 520, 521, 521, 521, 521, 521, 522, 522, 522, 522, 522, 523, 523, 523, 523, 523, 524, 524, 524, 524, 524, 525, 525, 525, 525, 525, 526, 526, 526, 526, 526, 527, 527, 527, 527, 527, 528, 528, 528, 528, 528, 529, 529, 529, 529, 529, 530, 530, 530, 530, 530, 531, 531, 531, 531, 531, 532, 532, 532, 532, 532, 533, 533, 533, 533, 533, 534, 534, 534, 534, 534, 535, 535, 535, 535, 535, 536, 536, 536, 536, 536, 537, 537, 537, 537, 537, 538, 538, 538, 538, 538, 539, 539, 539, 539, 539, 540, 540, 540, 540, 540, 541, 541, 541, 541, 541, 542, 542, 542, 542, 542, 543, 543, 543, 543, 543, 544, 544, 544, 544, 544, 545, 545, 545, 545, 545, 546, 546, 546, 546, 546, 547, 547, 547, 547, 547, 548, 548, 548, 548, 548, 549, 549, 549, 549, 549, 550, 550, 550, 550, 550, 551, 551, 551, 551, 551, 552, 552, 552, 552, 552, 553, 553, 553, 553, 553, 554, 554, 554, 554, 554, 555, 555, 555, 555, 555, 556, 556, 556, 556, 556, 557, 557, 557, 557, 557, 558, 558, 558, 558, 558, 559, 559, 559, 559, 559, 560, 560, 560, 560, 560, 561, 561, 561, 561, 561, 562, 562, 562, 562, 562, 563, 563, 563, 563, 563, 564, 564, 564, 564, 564, 565, 565, 565, 565, 565, 566, 566, 566, 566, 566, 567, 567, 567, 567, 567, 568, 568, 568, 568, 568, 569, 569, 569, 569, 569, 570, 570, 570, 570, 570, 571, 571, 571, 571, 571, 572, 572, 572, 572, 572, 573, 573, 573, 573, 573, 574, 574, 574, 574, 574, 575, 575, 575, 575, 575, 576, 576, 576, 576, 576, 577, 577, 577, 577, 577, 578, 578, 578, 578, 578, 579, 579, 579, 579, 579, 580, 580, 580, 580, 580, 581, 581, 581, 581, 581, 582, 582, 582, 582, 582, 583, 583, 583, 583, 583, 584, 584, 584, 584, 584, 585, 585, 585, 585, 585, 586, 586, 586, 586, 586, 587, 587, 587, 587, 587, 588, 588, 588, 588, 588, 589, 589, 589, 589, 589, 590, 590, 590, 590, 590, 591, 591, 591, 591, 591, 592, 592, 592, 592, 592, 593, 593, 593, 593, 593, 594, 594, 594, 594, 594, 595, 595, 595, 595, 595, 596, 596, 596, 596, 596, 597, 597, 597, 597, 597, 598, 598, 598, 598, 598, 599, 599, 599, 599, 599, 600, 600, 600, 600, 600, 601, 601, 601, 601, 601, 602, 602, 602, 602, 602, 603, 603, 603, 603, 603, 604, 604, 604, 604, 604, 605, 605, 605, 605, 605, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610, 610, 610, 610, 610, 611, 611, 611, 611, 611, 612, 612, 612, 612, 612, 613, 613, 613, 613, 613, 614, 614, 614, 614, 614, 615, 615, 615, 615, 615, 616, 616, 616, 616, 616, 617, 617, 617, 617, 617, 618, 618, 618, 618, 618, 619, 619, 619, 619, 619, 620, 620, 620, 620, 620, 621, 621, 621, 621, 621, 622, 622, 622, 622, 622, 623, 623, 623, 623, 623, 624, 624, 624, 624, 624, 625, 625, 625, 625, 625, 626, 626, 626, 626, 626, 627, 627, 627, 627, 627, 628, 628, 628, 628, 628, 629, 629, 629, 629, 629, 630, 630, 630, 630, 630, 631, 631, 631, 631, 631, 632, 632, 632, 632, 632, 633, 633, 633, 633, 633, 634, 634, 634, 634, 634, 635, 635, 635, 635, 635, 636, 636, 636, 636, 636, 637, 637, 637, 637, 637, 638, 638, 638, 638, 638, 639, 639, 639, 639, 639, 640, 640, 640, 640, 640, 641, 641, 641, 641, 641, 642, 642, 642, 642, 642, 643, 643, 643, 643, 643, 644, 644, 644, 644, 644, 645, 645, 645, 645, 645, 646, 646, 646, 646, 646, 647, 647, 647, 647, 647, 648, 648, 648, 648, 648, 649, 649, 649, 649, 649, 650, 650, 650, 650, 650, 651, 651, 651, 651, 651, 652, 652, 652, 652, 652, 653, 653, 653, 653, 653, 654, 654, 654, 654, 654, 655, 655, 655, 655, 655, 656, 656, 656, 656, 656, 657, 657, 657, 657, 657, 658, 658, 658, 658, 658, 659, 659, 659, 659, 659, 660, 660, 660, 660, 660, 661, 661, 661, 661, 661, 662, 662, 662, 662, 662, 663, 663, 663, 663, 663, 664, 664, 664, 664, 664, 665, 665, 665, 665, 665, 666, 666, 666, 666, 666, 667, 667, 667, 667, 667, 668, 668, 668, 668, 668, 669, 669, 669, 669, 669, 670, 670, 670, 670, 670, 671, 671, 671, 671, 671, 672, 672, 672, 672, 672, 673, 673, 673, 673, 673, 674, 674, 674, 674, 674, 675, 675, 675, 675, 675, 676, 676, 676, 676, 676, 677, 677, 677, 677, 677, 678, 678, 678, 678, 678, 679, 679, 679, 679, 679, 680, 680, 680, 680, 680, 681, 681, 681, 681, 681, 682, 682, 682, 682, 682, 683, 683, 683, 683, 683, 684, 684, 684, 684, 684, 685, 685, 685, 685, 685, 686, 686, 686, 686, 686, 687, 687, 687, 687, 687, 688, 688, 688, 688, 688, 689, 689, 689, 689, 689, 690, 690, 690, 690, 690, 691, 691, 691, 691, 691, 692, 692, 692, 692, 692, 693, 693, 693, 693, 693, 694, 694, 694, 694, 694, 695, 695, 695, 695, 695, 696, 696, 696, 696, 696, 697, 697, 697, 697, 697, 698, 698, 698, 698, 698, 699, 699, 699, 699, 699, 700, 700, 700, 700, 700, 701, 701, 701, 701, 701, 702, 702, 702, 702, 702, 703, 703, 703, 703, 703, 704, 704, 704, 704, 704, 705, 705, 705, 705, 705, 706, 706, 706, 706, 706, 707, 707, 707, 707, 707, 708, 708, 708, 708, 708, 709, 709, 709, 709, 709, 710, 710, 710, 710, 710, 711, 711, 711, 711, 711, 712, 712, 712, 712, 712, 713, 713, 713, 713, 713, 714, 714, 714, 714, 714, 715, 715, 715, 715, 715, 716, 716, 716, 716, 716, 717, 717, 717, 717, 717, 718, 718, 718, 718, 718, 719, 719, 719, 719, 719, 720, 720, 720, 720, 720, 721, 721, 721, 721, 721, 722, 722, 722, 722, 722, 723, 723, 723, 723, 723, 724, 724, 724, 724, 724, 725, 725, 725, 725, 725, 726, 726, 726, 726, 726, 727, 727, 727, 727, 727, 728, 728, 728, 728, 728, 729, 729, 729, 729, 729, 730, 730, 730, 730, 730, 731, 731, 731, 731, 731, 732, 732, 732, 732, 732, 733, 733, 733, 733, 733, 734, 734, 734, 734, 734, 735, 735, 735, 735, 735, 736, 736, 736, 736, 736, 737, 737, 737, 737, 737, 738, 738, 738, 738, 738, 739, 739, 739, 739, 739, 740, 740, 740, 740, 740, 741, 741, 741, 741, 741, 742, 742, 742, 742, 742, 743, 743, 743, 743, 743, 744, 744, 744, 744, 744, 745, 745, 745, 745, 745, 746, 746, 746, 746, 746, 747, 747, 747, 747, 747, 748, 748, 748, 748, 748, 749, 749, 749, 749, 749, 750, 750, 750, 750, 750, 751, 751, 751, 751, 751, 752, 752, 752, 752, 752, 753, 753, 753, 753, 753, 754, 754, 754, 754, 754, 755, 755, 755, 755, 755, 756, 756, 756, 756, 756, 757, 757, 757, 757, 757, 758, 758, 758, 758, 758, 759, 759, 759, 759, 759, 760, 760, 760, 760, 760, 761, 761, 761, 761, 761, 762, 762, 762, 762, 762, 763, 763, 763, 763, 763, 764, 764, 764, 764, 764, 765, 765, 765, 765, 765, 766, 766, 766, 766, 766, 767, 767, 767, 767, 767, 768, 768, 768, 768, 768, 769, 769, 769, 769, 769, 770, 770, 770, 770, 770, 771, 771, 771, 771, 771, 772, 772, 772, 772, 772, 773, 773, 773, 773, 773, 774, 774, 774, 774, 774, 775, 775, 775, 775, 775, 776, 776, 776, 776, 776, 777, 777, 777, 777, 777, 778, 778, 778, 778, 778, 779, 779, 779, 779, 779, 780, 780, 780, 780, 780, 781, 781, 781, 781, 781, 782, 782, 782, 782, 782, 783, 783, 783, 783, 783, 784, 784, 784, 784, 784, 785, 785, 785, 785, 785, 786, 786, 786, 786, 786, 787, 787, 787, 787, 787, 788, 788, 788, 788, 788, 789, 789, 789, 789, 789, 790, 790, 790, 790, 790, 791, 791, 791, 791, 791, 792, 792, 792, 792, 792, 793, 793, 793, 793, 793, 794, 794, 794, 794, 794, 795, 795, 795, 795, 795, 796, 796, 796, 796, 796, 797, 797, 797, 797, 797, 798, 798, 798, 798, 798, 799, 799, 799, 799, 799, 800, 800, 800, 800, 800, 801, 801, 801, 801, 801, 802, 802, 802, 802, 802, 803, 803, 803, 803, 803, 804, 804, 804, 804, 804, 805, 805, 805, 805, 805, 806, 806, 806, 806, 806, 807, 807, 807, 807, 807, 808, 808, 808, 808, 808, 809, 809, 809, 809, 809, 810, 810, 810, 810, 810, 811, 811, 811, 811, 811, 812, 812, 812, 812, 812, 813, 813, 813, 813, 813, 814, 814, 814, 814, 814, 815, 815, 815, 815, 815, 816, 816, 816, 816, 816, 817, 817, 817, 817, 817, 818, 818, 818, 818, 818, 819, 819, 819, 819, 819, 820, 820, 820, 820, 820, 821, 821, 821, 821, 821, 822, 822, 822, 822, 822, 823, 823, 823, 823, 823, 824, 824, 824, 824, 824, 825, 825, 825, 825, 825, 826, 826, 826, 826, 826, 827, 827, 827, 827, 827, 828, 828, 828, 828, 828, 829, 829, 829, 829, 829, 830, 830, 830, 830, 830, 831, 831, 831, 831, 831, 832, 832, 832, 832, 832, 833, 833, 833, 833, 833, 834, 834, 834, 834, 834, 835, 835, 835, 835, 835, 836, 836, 836, 836, 836, 837, 837, 837, 837, 837, 838, 838, 838, 838, 838, 839, 839, 839, 839, 839, 840, 840, 840, 840, 840, 841, 841, 841, 841, 841, 842, 842, 842, 842, 842, 843, 843, 843, 843, 843, 844, 844, 844, 844, 844, 845, 845, 845, 845, 845, 846, 846, 846, 846, 846, 847, 847, 847, 847, 847, 848, 848, 848, 848, 848, 849, 849, 849, 849, 849, 850, 850, 850, 850, 850, 851, 851, 851, 851, 851, 852, 852, 852, 852, 852, 853, 853, 853, 853, 853, 854, 854, 854, 854, 854, 855, 855, 855, 855, 855, 856, 856, 856, 856, 856, 857, 857, 857, 857, 857, 858, 858, 858, 858, 858, 859, 859, 859, 859, 859, 860, 860, 860, 860, 860, 861, 861, 861, 861, 861, 862, 862, 862, 862, 862, 863, 863, 863, 863, 863, 864, 864, 864, 864, 864, 865, 865, 865, 865, 865, 866, 866, 866, 866, 866, 867, 867, 867, 867, 867, 868, 868, 868, 868, 868, 869, 869, 869, 869, 869, 870, 870, 870, 870, 870, 871, 871, 871, 871, 871, 872, 872, 872, 872, 872, 873, 873, 873, 873, 873, 874, 874, 874, 874, 874, 875, 875, 875, 875, 875, 876, 876, 876, 876, 876, 877, 877, 877, 877, 877, 878, 878, 878, 878, 878, 879, 879, 879, 879, 879, 880, 880, 880, 880, 880, 881, 881, 881, 881, 881, 882, 882, 882, 882, 882, 883, 883, 883, 883, 883, 884, 884, 884, 884, 884, 885, 885, 885, 885, 885, 886, 886, 886, 886, 886, 887, 887, 887, 887, 887, 888, 888, 888, 888, 888, 889, 889, 889, 889, 889, 890, 890, 890, 890, 890, 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 893, 893, 893, 893, 894, 894, 894, 894, 894, 895, 895, 895, 895, 895, 896, 896, 896, 896, 896, 897, 897, 897, 897, 897, 898, 898, 898, 898, 898, 899, 899, 899, 899, 899, 900, 900, 900, 900, 900, 901, 901, 901, 901, 901, 902, 902, 902, 902, 902, 903, 903, 903, 903, 903, 904, 904, 904, 904, 904, 905, 905, 905, 905, 905, 906, 906, 906, 906, 906, 907, 907, 907, 907, 907, 908, 908, 908, 908, 908, 909, 909, 909, 909, 909, 910, 910, 910, 910, 910, 911, 911, 911, 911, 911, 912, 912, 912, 912, 912, 913, 913, 913, 913, 913, 914, 914, 914, 914, 914, 915, 915, 915, 915, 915, 916, 916, 916, 916, 916, 917, 917, 917, 917, 917, 918, 918, 918, 918, 918, 919, 919, 919, 919, 919, 920, 920, 920, 920, 920, 921, 921, 921, 921, 921, 922, 922, 922, 922, 922, 923, 923, 923, 923, 923, 924, 924, 924, 924, 924, 925, 925, 925, 925, 925, 926, 926, 926, 926, 926, 927, 927, 927, 927, 927, 928, 928, 928, 928, 928, 929, 929, 929, 929, 929, 930, 930, 930, 930, 930, 931, 931, 931, 931, 931, 932, 932, 932, 932, 932, 933, 933, 933, 933, 933, 934, 934, 934, 934, 934, 935, 935, 935, 935, 935, 936, 936, 936, 936, 936, 937, 937, 937, 937, 937, 938, 938, 938, 938, 938, 939, 939, 939, 939, 939, 940, 940, 940, 940, 940, 941, 941, 941, 941, 941, 942, 942, 942, 942, 942, 943, 943, 943, 943, 943, 944, 944, 944, 944, 944, 945, 945, 945, 945, 945, 946, 946, 946, 946, 946, 947, 947, 947, 947, 947, 948, 948, 948, 948, 948, 949, 949, 949, 949, 949, 950, 950, 950, 950, 950, 951, 951, 951, 951, 951, 952, 952, 952, 952, 952, 953, 953, 953, 953, 953, 954, 954, 954, 954, 954, 955, 955, 955, 955, 955, 956, 956, 956, 956, 956, 957, 957, 957, 957, 957, 958, 958, 958, 958, 958, 959, 959, 959, 959, 959, 960, 960, 960, 960, 960, 961, 961, 961, 961, 961, 962, 962, 962, 962, 962, 963, 963, 963, 963, 963, 964, 964, 964, 964, 964, 965, 965, 965, 965, 965, 966, 966, 966, 966, 966, 967, 967, 967, 967, 967, 968, 968, 968, 968, 968, 969, 969, 969, 969, 969, 970, 970, 970, 970, 970, 971, 971, 971, 971, 971, 972, 972, 972, 972, 972, 973, 973, 973, 973, 973, 974, 974, 974, 974, 974, 975, 975, 975, 975, 975, 976, 976, 976, 976, 976, 977, 977, 977, 977, 977, 978, 978, 978, 978, 978, 979, 979, 979, 979, 979, 980, 980, 980, 980, 980, 981, 981, 981, 981, 981, 982, 982, 982, 982, 982, 983, 983, 983, 983, 983, 984, 984, 984, 984, 984, 985, 985, 985, 985, 985, 986, 986, 986, 986, 986, 987, 987, 987, 987, 987, 988, 988, 988, 988, 988, 989, 989, 989, 989, 989, 990, 990, 990, 990, 990, 991, 991, 991, 991, 991, 992, 992, 992, 992, 992, 993, 993, 993, 993, 993, 994, 994, 994, 994, 994, 995, 995, 995, 995, 995, 996, 996, 996, 996, 996, 997, 997, 997, 997, 997, 998, 998, 998, 998, 998, 999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1001, 1001, 1001, 1001, 1001, 1002, 1002, 1002, 1002, 1002, 1003, 1003, 1003, 1003, 1003, 1004, 1004, 1004, 1004, 1004, 1005, 1005, 1005, 1005, 1005, 1006, 1006, 1006, 1006, 1006, 1007, 1007, 1007, 1007, 1007, 1008, 1008, 1008, 1008, 1008, 1009, 1009, 1009, 1009, 1009, 1010, 1010, 1010, 1010, 1010, 1011, 1011, 1011, 1011, 1011, 1012, 1012, 1012, 1012, 1012, 1013, 1013, 1013, 1013, 1013, 1014, 1014, 1014, 1014, 1014, 1015, 1015, 1015, 1015, 1015, 1016, 1016, 1016, 1016, 1016, 1017, 1017, 1017, 1017, 1017, 1018, 1018, 1018, 1018, 1018, 1019, 1019, 1019, 1019, 1019, 1020, 1020, 1020, 1020, 1020, 1021, 1021, 1021, 1021, 1021, 1022, 1022, 1022, 1022, 1022, 1023, 1023, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1024], [4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0  …  -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0], 1024, 1024), [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], [0.0 0.03842943919353914 … 0.15224093497742697 0.03842943919353936; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785; … ; 0.15224093497742697 0.1906703741709661 … 0.30448186995485393 0.19067037417096633; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785])

One of the benefits of the Bayesian approach is that we can fit for the hyperparameters of our prior/regularizers unlike traditional RML appraoches. To construct this heirarchical prior we will first make a map that takes in our regularizer hyperparameters and returns the image prior given those hyperparameters.

fmap = let crcache=crcache
    x->GaussMarkovRandomField(x, 1.0, crcache)
end
#1 (generic function with 1 method)

Now we can finally form our image prior. For this we use a heirarchical prior where the inverse correlation length is given by a Half-Normal distribution whose peak is at zero and standard deviation is 0.1/rat where recall rat is the beam size per pixel. For the variance of the random field we use another half normal prior with standard deviation 0.1. The reason we use the half-normal priors is to prefer "simple" structures. Gaussian Markov random fields are extremly flexible models, and to prevent overfitting it is common to use priors that penalize complexity. Therefore, we want to use priors that enforce similarity to our mean image. If the data wants more complexity then it will drive us away from the prior.

cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.01*rat)))
VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)

We can now form our model parameter priors. Like our other imaging examples, we use a Dirichlet prior for our image pixels. For the log gain amplitudes, we use the CalPrior which automatically constructs the prior for the given jones cache gcache.

prior = NamedDist(
         c = cprior,
         fg = Uniform(0.0, 1.0),
         σimg = truncated(Normal(0.0, 0.1); lower = 0.0),
         lgamp = CalPrior(distamp, gcache),
         gphase = CalPrior(distphase, gcachep),
        )
VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)
, Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0), CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)

jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778  …  0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)

jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)

Putting it all together we form our likelihood and posterior objects for optimization and sampling.

lklhd = RadioLikelihood(sky, instrument, dvis; skymeta=metadata, instrumentmeta=metadata)
post = Posterior(lklhd, prior)
Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Comrade.ModelMetadata{typeof(Main.instrument), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#26#27"{Float64, Vector{Float64}}, Vector{ComplexF64}}}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}}(RadioLikelihood
	Number of data products: 1
, VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)
, Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0), CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)

jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778  …  0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)

jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)
)

Reconstructing the Image and Instrument Effects

To sample from this posterior, it is convenient to move from our constrained parameter space to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This is done using the asflat function.

tpost = asflat(post)
ndim = dimension(tpost)
1364

Our Posterior and TransformedPosterior objects satisfy the LogDensityProblems interface. This allows us to easily switch between different AD backends and many of Julia's statistical inference packages use this interface as well.

using LogDensityProblemsAD
using Zygote
gtpost = ADgradient(Val(:Zygote), tpost)
x0 = randn(rng, ndim)
LogDensityProblemsAD.logdensity_and_gradient(gtpost, x0)
(-1.0889677227706836e8, [3.9714863653464376, -41.829408532001366, 69.59557261663501, 4.569805366219422, 10.789295400079082, 240.91028613995235, 39.14697994202466, 42.2142944805375, 50.171331063617025, 69.99422943458495  …  3666.7681059308334, 11816.75911472515, 1587.1490272240208, -343.3074914007036, 6760.359870612455, 349.8400750594293, 17534.646730431836, -7352.639076548009, -181.3718982383252, 410.63933095430605])

We can now also find the dimension of our posterior or the number of parameters we are going to sample.

Warning

This can often be different from what you would expect. This is especially true when using angular variables where we often artificially increase the dimension of the parameter space to make sampling easier.

To initialize our sampler we will use optimize using LBFGS

using ComradeOptimization
using OptimizationOptimJL
f = OptimizationFunction(tpost, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing)
ℓ = logdensityof(tpost)
sol = solve(prob, LBFGS(), maxiters=4_000, g_tol=1e-1);

Now transform back to parameter space

xopt = transform(tpost, sol.u)
(c = (params = [-6.037456250054449e-6 -6.571354045485335e-6 … -1.85679528330551e-6 -2.291977091024625e-6; -1.0648811374208013e-5 -1.7981714076150314e-5 … -1.1167696480467717e-5 -6.0748012291175495e-6; … ; 9.157073932541369e-6 4.918747579945945e-6 … 6.55279425861554e-6 2.0301535450304157e-6; 2.956326326988618e-6 1.107977057503808e-6 … 6.472800399555038e-6 1.689377636932908e-6], hyperparams = 0.056195917416529226), fg = 0.003549803518697557, σimg = 3.1071521396538953, lgamp = [-0.009113077870190382, 0.009107168806961666, -0.8878392223496306, 0.13651961509739227, 0.014518947611761137, -0.014031849197920259, -0.316678910079671, 0.24192237915862175, 0.005584388045765145, -0.005197415859705387  …  -0.29303878993604465, -0.0006552411035840428, -1.0959177220925187, -0.003013688240126423, 0.026843947617235218, -0.026107816599595716, -0.2822707440520687, 0.000568849136452204, -1.1550018292416246, -0.004406101066821852], gphase = [-0.8142655675167031, -3.1364555519407404, 2.351460373182827, -0.8711498011840125, -3.038854928403769, 2.157523854335552, -0.9285480596058518, -2.952244171814421, 1.9680821210601833, -0.9773334792754164  …  -0.5846216803138383, -0.4926521978100723, -0.5149517295308924, -0.6901891406444994, 1.2564547496971952, -0.5238549827932711, -0.6053695673999937, -0.48455476766515043, -0.9024513472709461, 1.3044807696191534])
Warning

Fitting gains tends to be very difficult, meaning that optimization can take a lot longer. The upside is that we usually get nicer images.

First we will evaluate our fit by plotting the residuals

using Plots
residual(vlbimodel(post, xopt), dvis)
Example block output

These look reasonable, although there may be some minor overfitting. This could be improved in a few ways, but that is beyond the goal of this quick tutorial. Plotting the image, we see that we have a much cleaner version of the closure-only image from Imaging a Black Hole using only Closure Quantities.

import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), fovx, fovy, 128, 128)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot)

Because we also fit the instrument model, we can inspect their parameters. To do this, Comrade provides a caltable function that converts the flattened gain parameters to a tabular format based on the time and its segmentation.

gt = Comrade.caltable(gcachep, xopt.gphase)
plot(gt, layout=(3,3), size=(600,500))
Example block output

The gain phases are pretty random, although much of this is due to us picking a random reference station for each scan.

Moving onto the gain amplitudes, we see that most of the gain variation is within 10% as expected except LMT, which has massive variations.

gt = Comrade.caltable(gcache, exp.(xopt.lgamp))
plot(gt, layout=(3,3), size=(600,500))
Example block output

To sample from the posterior, we will use HMC, specifically the NUTS algorithm. For information about NUTS, see Michael Betancourt's notes.

Note

For our metric, we use a diagonal matrix due to easier tuning

However, due to the need to sample a large number of gain parameters, constructing the posterior is rather time-consuming. Therefore, for this tutorial, we will only do a quick preliminary run, and any posterior inferences should be appropriately skeptical.

using ComradeAHMC
metric = DiagEuclideanMetric(ndim)
chain, stats = sample(rng, post, AHMC(;metric, autodiff=Val(:Zygote)), 700; nadapts=500, init_params=xopt)
(NamedTuple{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{NamedTuple{(:params, :hyperparams), Tuple{Matrix{Float64}, Float64}}, Float64, Float64, Vector{Float64}, Vector{Float64}}}[(c = (params = [-0.001576473654468533 0.007447526608483789 … 0.01170705016326149 -0.007005803664701813; 0.007163398298646599 0.013367675949110057 … 0.0001579602080998415 -0.0033051406764038264; … ; -0.005643781438394142 0.004440087675231427 … -0.004533397876913021 0.0016239710571917; -0.009682768008299768 0.010530390614365354 … -0.0035480079606065758 -0.003306247524249095], hyperparams = 0.0562680467523908), fg = 0.0035896791710454566, σimg = 3.1212668715336096, lgamp = [-0.009838565836800562, 0.009245813397830637, -0.8916877002951414, 0.13431256103024947, 0.02114880609643208, -0.02104395619614467, -0.32493354003428776, 0.2292686122979385, 0.004403670989606053, -0.003916563038231321  …  -0.2943496878027941, 0.0019100561178617362, -1.1030655592825191, 0.0016072678897755678, 0.02440463945920735, -0.024019566237323587, -0.28426035726311377, -0.006905917741296118, -1.1589013405299586, 0.003118249639101413], gphase = [-0.8173284134425431, -3.120496658804397, 2.3597920296228074, -0.8685770594169238, -3.050684457990954, 2.1663650203931493, -0.9280517378643218, -2.953797960904017, 1.9714725593121756, -0.978753097320341  …  -0.5869695513173957, -0.4992328258229167, -0.5013654966677868, -0.6904707013965257, 1.26744067387979, -0.5219069075973074, -0.6030348852808273, -0.4862872214642622, -0.9050541659918923, 1.2985782475644916]), (c = (params = [-0.001576473654468533 0.007447526608483789 … 0.01170705016326149 -0.007005803664701813; 0.007163398298646599 0.013367675949110057 … 0.0001579602080998415 -0.0033051406764038264; … ; -0.005643781438394142 0.004440087675231427 … -0.004533397876913021 0.0016239710571917; -0.009682768008299768 0.010530390614365354 … -0.0035480079606065758 -0.003306247524249095], hyperparams = 0.0562680467523908), fg = 0.0035896791710454566, σimg = 3.1212668715336096, lgamp = [-0.009838565836800562, 0.009245813397830637, -0.8916877002951414, 0.13431256103024947, 0.02114880609643208, -0.02104395619614467, -0.32493354003428776, 0.2292686122979385, 0.004403670989606053, -0.003916563038231321  …  -0.2943496878027941, 0.0019100561178617362, -1.1030655592825191, 0.0016072678897755678, 0.02440463945920735, -0.024019566237323587, -0.28426035726311377, -0.006905917741296118, -1.1589013405299586, 0.003118249639101413], gphase = [-0.8173284134425431, -3.120496658804397, 2.3597920296228074, -0.8685770594169238, -3.050684457990954, 2.1663650203931493, -0.9280517378643218, -2.953797960904017, 1.9714725593121756, -0.978753097320341  …  -0.5869695513173957, -0.4992328258229167, -0.5013654966677868, -0.6904707013965257, 1.26744067387979, -0.5219069075973074, -0.6030348852808273, -0.4862872214642622, -0.9050541659918923, 1.2985782475644916]), (c = (params = [-0.11446803346348784 -0.013823627054498609 … -0.17011746213487497 -0.12919708408893632; 0.012498207367954888 -0.024006053934599026 … -0.049954591436122006 -0.023953044628096803; … ; -0.08883473867667067 0.0010521703398028758 … -0.10507474900063955 -0.09342913521951575; -0.10905941867606325 -0.035639270894171134 … -0.058276784298745575 0.06184044650698433], hyperparams = 0.1192679481314917), fg = 0.0028764812445817222, σimg = 2.545773460104071, lgamp = [0.0020960537497268714, -0.005059934546357328, -0.8852541327152603, 0.11104960516010773, -0.0015784983278665206, 0.003418795959386172, -0.3114903300449294, 0.24856687972621605, 0.00016154551654989696, -0.001005216391175049  …  -0.30045287760480827, -0.0065261992320247856, -1.1148678270518513, -0.011258760341383005, 0.026127588148970166, -0.024961911849735348, -0.27699194783655307, -0.0007544533728501114, -1.165545760614164, -0.00324373035909087], gphase = [-0.8187041767986699, 3.128401096173404, 2.3799980223111774, -0.8665248865498782, -3.0356059419283463, 2.1853183879209492, -0.9293000283621785, -2.948883082533599, 1.9961476088431325, -0.9776084088361938  …  -0.5811378498738728, -0.5001188375581221, -0.5366041012227611, -0.6918023176995408, 1.2348174902075597, -0.5211135553227288, -0.6108145559146654, -0.5176133441381178, -0.8965181620200161, 1.2627068249931157]), (c = (params = [0.04240657902355421 0.09721626271064752 … -0.30580927430051746 0.0015614730422718738; -0.025444717984635355 -0.19073828346203248 … -0.16675721272805255 0.04997049296390746; … ; 0.032023734250238835 -0.040643711248192514 … -0.06568473075331102 -0.10972421333165011; -0.14577603288446217 -0.09083404430866811 … 0.010700425757568978 0.11684502175445267], hyperparams = 0.16969685946250407), fg = 0.0035322909355134506, σimg = 2.021740793112091, lgamp = [0.0016970215064956668, -0.0018016220142085437, -0.9899929300068033, 0.07177506621966387, -0.005326039383352655, 0.0013673405163830872, -0.3857619361018142, 0.20609949631765767, -0.008944835920071016, 0.002127326341638535  …  -0.2915070534717762, -0.017669473947270316, -1.1144170997837006, 0.0003050684323877977, 0.00925161951600908, -0.013309871382162776, -0.26952029814762013, 0.008427956763776258, -1.1689893826428135, -0.01059757185348734], gphase = [-0.8151832977618284, -3.08732286523961, 2.3713433292846924, -0.8671237132000987, -3.0040271475541758, 2.182338175008595, -0.9333964588391551, -2.9432204991252213, 1.9949612971681223, -0.9770670252862125  …  -0.5871061683387708, -0.4818633020904415, -0.5786657258967866, -0.6848039187476445, 1.1934832407854574, -0.52741049893656, -0.6081019482761871, -0.5202401236986663, -0.896948569127936, 1.262682832585686]), (c = (params = [-0.12013742748287824 -0.13365567639338116 … 0.3598059770164582 0.05017767860751387; 0.040366494391990006 0.22576426001985903 … 0.15507351167638653 -0.11609004098464108; … ; 0.026889458351487636 -0.007491090660486707 … 0.015655333249512224 0.16267485944478055; 0.23072674990114622 0.07558070500025997 … -0.07881021506696452 -0.08247667366084926], hyperparams = 0.24784509057153403), fg = 0.011663229266328314, σimg = 1.5779106960795337, lgamp = [0.008509729806815974, -0.009728758243039265, -1.106285809896934, -0.04577871893750097, -0.004073507646833414, 0.004050929113822096, -0.5032037922324583, 0.10960570818431027, -0.018906502891660945, 0.019335878643881595  …  -0.2701537214214951, 0.010953330363942764, -1.1398579170298822, -0.014007370842038814, -0.007711652054149772, 0.004509078311824885, -0.25239484783574756, -0.0057824949771027555, -1.2085878680350415, 0.003913987061647158], gphase = [-0.8129493282072382, -3.131862404963679, 2.2312082137624003, -0.8693328006809525, -3.018655598729151, 2.031052052309857, -0.9287579997191622, -2.9546164523093985, 1.8411936649548744, -0.9738069993096858  …  -0.5814059428834325, -0.49028889570948814, -0.4141706958038609, -0.6841441340396814, 1.3527765214420708, -0.5307066259419684, -0.6115618379687473, -0.3556253811431712, -0.895893012005879, 1.4316868178024897]), (c = (params = [0.042526861146469835 0.14179683550190247 … -0.4360436652854706 -0.1281617752720419; -0.04241126228516675 -0.35546858602237813 … -0.18821112929945338 0.06360504503705684; … ; -0.05373973892563489 -0.013650292966841147 … -0.043368663530207095 -0.20556537588975424; -0.2654919264437389 -0.1156053158282567 … 0.08400147374524085 0.06438839636875016], hyperparams = 0.24614521745356177), fg = 0.022413764714074805, σimg = 1.435106806323219, lgamp = [-0.03494585823975815, 0.03322370968389133, -1.0556509523248134, 0.011991828629396982, -0.008569837044652587, 0.01285689054978069, -0.49957301651013103, 0.10985602448899186, -0.02579394414411876, 0.027625891838015487  …  -0.2649204023971786, -0.006719381924882974, -1.1388387692805282, -0.009882117090429157, -0.01698875694040232, 0.011746043056718714, -0.23092891941714097, -0.009096005451401257, -1.1908121513969059, 0.007109213668632348], gphase = [-0.8089437700135729, -3.022024744462354, 2.33622097821036, -0.8727379442352616, -2.9978448073247037, 2.1463388528893934, -0.9306969224437515, -2.9243665349366936, 1.9689494065642388, -0.9769412509971718  …  -0.5826210029211071, -0.48299309217029823, -0.4752059034311415, -0.6938386197756496, 1.291791053757474, -0.5172805150912648, -0.6083335664748106, -0.460850506360218, -0.9021107029712189, 1.3318573330833148]), (c = (params = [0.03667764323982257 -0.2888708981335546 … -0.14581370129050775 0.3607943690537443; -0.2545624779999389 -0.5757975260916753 … -0.29015524742285353 0.05053301680109258; … ; 0.13203343843282714 -0.471461665027403 … -0.15716934032599736 -0.11738360567706052; 0.010599267716610052 -0.3180302124581771 … 0.024349135732820052 0.5875459499588289], hyperparams = 0.3354562296532181), fg = 0.02272286140420574, σimg = 1.6317535007302697, lgamp = [-0.016333368839066366, 0.01776176002466547, -1.1151325014309443, 0.04072815854447249, -0.046962775696581735, 0.0460300707006579, -0.4640433684076863, 0.20187520260852831, -0.021765887850908494, 0.023178209561426088  …  -0.2457090781468747, -0.006660547251218679, -1.129234105835501, 0.007731822363585434, -0.011204732964525517, 0.00788427151309323, -0.26683323356497746, -0.015790124977789113, -1.18964710644862, 0.007249528374559371], gphase = [-0.8249294695268473, -3.021823986052246, 2.2892907735137573, -0.872314865088574, -2.9812849979022666, 2.0866589224505523, -0.9258838594791483, -2.9160763831710454, 1.9047690978888407, -0.9779302543896358  …  -0.5823638960670379, -0.49914233800661806, -0.42597395948454525, -0.678261728447934, 1.3370413100260317, -0.5229046108635123, -0.6003185831738775, -0.3849053959155854, -0.8632728725711293, 1.394999840317212]), (c = (params = [0.39249991277975427 -0.18929678988337018 … 0.6002278106273515 -0.19946073116919155; 0.17439396492271014 0.2048036119203854 … 0.2849708665868768 -0.07327518945667141; … ; -0.782703849996918 0.13313367537777432 … -0.09446382039423462 -0.11974934982896814; -0.2821450819581035 -0.06050735876230002 … -0.36534310472035686 -0.6341208296811758], hyperparams = 0.5229573765166583), fg = 0.032395989741656014, σimg = 1.2982276361552199, lgamp = [0.020720845145585857, -0.020993333705497212, -1.2123624819790801, -0.038251145535322294, -0.05590095314446232, 0.05329713668233851, -0.4899897498326019, 0.15930088334235012, -0.010259603186801157, 0.011420411294855598  …  -0.24729059060115297, -0.011596530088204382, -1.1261221234533356, 0.0021482387306447884, 0.02079422704609887, -0.024342217015182462, -0.26181156226052515, -0.02017095093970658, -1.1997064130704247, 0.007685786660610443], gphase = [-0.8186756114973744, -3.0560690679770683, 2.0923751401946733, -0.8721887928075787, -2.941925109134828, 1.9014624863748932, -0.9271515553743757, -2.86694908554747, 1.7199430427301803, -0.9748572516527815  …  -0.5844454733624992, -0.4630257189749749, -0.2778160043518557, -0.6568272977149976, 1.4976348129040538, -0.5249465426738418, -0.5553722107375048, -0.26018240183222274, -0.8759454263688904, 1.5405977933270603]), (c = (params = [-0.2403766436175168 0.3742725247358248 … 0.9100454466570612 0.7813406659869196; -0.1897401860149622 -0.019707363108415518 … -0.47284612447202806 0.33620156048160926; … ; -0.4741753234819241 0.1271896730488913 … 0.6049395868839416 0.10357770620939273; 0.05966007548306427 0.2521557550447108 … 0.1549967364453539 0.8035698038999886], hyperparams = 0.5610867561023429), fg = 0.03367628269714859, σimg = 1.3248845049340854, lgamp = [-0.004392527160900283, 0.00442075729445633, -1.1232532316292791, -0.030986576058934377, -0.046569088923891645, 0.048969364860778374, -0.5011506851799414, 0.14230599956023804, -0.010290129099326785, 0.007046360790064821  …  -0.23409046854370127, -0.0070111811280444795, -1.1043944128299414, -0.0068020045710442, -0.00596599566763031, 0.006408931039551745, -0.2531243808489334, 0.009647940970052733, -1.1661575818226775, -0.019480960310184356], gphase = [-0.8088563675171683, -3.011432236525387, 2.1878341020042593, -0.8707722848236895, -2.9531903686382432, 2.0127187815418006, -0.9263443061164652, -2.8731760551745507, 1.8289354610539221, -0.9767566701613045  …  -0.5853471586505526, -0.47162055368861083, -0.35242882560738364, -0.6378276052279811, 1.4158413972787403, -0.5233092984508445, -0.5759967168999022, -0.3137066614360691, -0.8355862486658542, 1.4684759762002104]), (c = (params = [0.3329379869827929 -0.48397675960945163 … -0.7693848677271717 -0.4692802274405943; -0.09499609026411664 -0.5022692992218208 … 0.9811535227535839 -0.061322802435992134; … ; 0.513106548860822 -0.03576717161812522 … -0.8241125274983153 -0.33326483055744516; 0.06000264293029741 -0.0013314048610602123 … -0.11596090898647303 -0.8580096292658689], hyperparams = 1.1289945929849852), fg = 0.47830150254455917, σimg = 1.2872369327201647, lgamp = [-0.03702466336382978, 0.03414583484589659, -0.7930177722233172, 0.043864489377270485, 0.004006637090615503, -0.0005450994958012636, -0.24919352001423473, 0.13236971219943322, -0.019110006790296508, 0.014394055006147844  …  0.004112142906542318, 0.01378713885581914, -0.8182383112337027, -0.008180127658734777, 0.021374659010452865, -0.01928930548669743, -0.006860249652773929, 0.0049671709471819005, -0.887295956968624, -0.005486691820309631], gphase = [-0.8186779095392948, -2.8674112583049376, 2.3461259742555454, -0.8729212501097275, -2.821888583516067, 2.169061826646109, -0.9268952047402386, -2.7823734175848367, 1.9829420580854553, -0.9800234988982228  …  -0.584105424929569, -0.3230890441286844, -0.26729953907379805, -0.5681247351144338, 1.4992893833104697, -0.5216425530023105, -0.4644351848201822, -0.2177115739773855, -0.752737897772711, 1.5773975274900174])  …  (c = (params = [-0.22729630457830385 0.7177195610019145 … 0.23117266935303155 0.8287702845681999; -0.03922959153730689 -0.5469540355902806 … -0.08899422503163644 0.05774367837124978; … ; 1.4416508818751574 0.784880832249737 … 2.4362599451131213 1.4192332227126312; -0.8476170125997223 0.3494835153637129 … 0.7989053674420142 0.5162414913154592], hyperparams = 13.453515967535315), fg = 0.3481459779422762, σimg = 0.9276577426536813, lgamp = [-0.034095013838816274, 0.03441531794385418, -0.9559624809926947, 0.030576394080644307, -0.02449347264417922, 0.02003559159636763, -0.35745887262651155, 0.195879250192966, -0.052008778673326886, 0.051816381777267495  …  -0.07990941751257628, 0.008568216804399797, -0.9348173046626438, -0.01387287549435151, 0.015787002994990324, -0.017311358943677562, -0.09050678924282721, 0.002060003611809816, -0.9746379743570998, -0.0016616786671206196], gphase = [-0.8126856333995779, -2.952867551583251, 0.7556613473859807, -0.8653983032696606, -2.8201872375859143, 0.5509236783664712, -0.9242171565582951, -2.749487015280738, 0.3726540003365169, -0.9827518704601854  …  -0.5813415524117573, -0.26367007552445476, 0.9818430048194294, -0.5094390625504188, 2.746085923321215, -0.5259560974381549, -0.4238714843146229, 0.9891421908019208, -0.733858877272575, 2.781791237076261]), (c = (params = [-0.16264841382145012 0.19229053244857996 … 0.5678204750439416 -0.23259678278748644; -0.6113982049836926 -0.7501178497668701 … 0.6229043681522815 -0.31635514934728937; … ; 0.6161195113352999 1.1722726011760203 … -0.031931903550287896 -0.02885756630899227; 1.5583781731767652 0.1119520860394075 … 0.4239366286453495 0.5226200982385968], hyperparams = 29.596705725481936), fg = 0.32037537779802994, σimg = 1.0489729105010208, lgamp = [-0.04203795231425661, 0.0447934584998382, -0.9834133358936458, 0.002598782982534062, -0.016453881861316882, 0.013576403108344189, -0.3685291382619457, 0.12822538540928857, -0.0033527267089412897, 0.0002666903096856762  …  -0.08679853363906942, -0.00928624261045769, -0.9414358838764266, 0.004334817230557164, 0.012320933167741922, -0.011331729892989582, -0.12458589586726468, -0.0013065721369138706, -1.0380612108587395, 0.0024702855909443664], gphase = [-0.8168399260924719, -2.9670171736680793, 1.0103740156856564, -0.8668616245300824, -2.887748715978341, 0.7870810265694324, -0.9328566934191805, -2.8152643524493444, 0.5845125934552582, -0.9841702324247247  …  -0.5889537491248853, -0.28883201086996413, 1.0555361375822636, -0.5576293770314716, 2.8290163329826905, -0.5277245789016936, -0.4402424039338409, 1.112668346882871, -0.7898649099561815, 2.8848120743372423]), (c = (params = [0.9556427666411855 0.9030028050930629 … 0.7324272012425317 1.2488855641091678; 1.620447538100932 0.6703238783949531 … 0.6083875659805043 1.1143916959525908; … ; 0.6433845734228876 0.999802611185972 … 1.5448965913481525 1.9174734402314126; -0.4158995559738227 0.5491556513157886 … 0.7899038506080188 0.20008069710798987], hyperparams = 15.20376735831896), fg = 0.307401032754837, σimg = 0.9267717378361005, lgamp = [-0.05559229428717584, 0.05291645904476581, -0.9568989927363047, -0.0076269047482570625, -0.008063026364300977, 0.0026261305703679377, -0.3698264940309769, 0.12562761985716966, -0.04258198812760945, 0.04079483027893814  …  -0.11416847522582146, 0.010288378684580489, -0.9378480670984819, -0.02322944863012248, -0.007514963791054445, 0.002656254882386254, -0.11005483373072941, 0.012384699765794425, -0.9695194689818627, -0.009637636110720648], gphase = [-0.811318066277044, -3.079077272488328, 0.9015037417273272, -0.8731358177884142, -2.988366547249762, 0.6946377447197728, -0.9327399445129426, -2.899450587245803, 0.5007102848602433, -0.9799781720448196  …  -0.5869747640621799, -0.3095454837426049, 0.9793904998722988, -0.563694131820663, 2.7539926017036365, -0.5180432561910668, -0.44071486709884267, 1.002434713084724, -0.7868324329433679, 2.775813409929233]), (c = (params = [0.7483403439603289 0.4063989538421792 … 1.1240449464445936 0.8818716827988571; -0.5222559195074283 -0.15629877499006445 … 1.3869186440042978 0.3167875423496312; … ; 1.1530194388416726 0.5807994961888951 … -0.20039197815752488 -0.07058023983943149; 2.491868549976821 0.5988474990689867 … 0.8869423752324982 1.5221340320797234], hyperparams = 27.952680362600347), fg = 0.3290893237663129, σimg = 0.956541307414417, lgamp = [-0.0013276440801088347, 0.008974500148812101, -1.0257242415087064, 0.029365392951845502, -0.0011151073705694282, 0.0011590114337681795, -0.4085083707946752, 0.16218397019219427, -0.015226219484989468, 0.016366685085706975  …  -0.06842245997339838, 0.017351249524479542, -0.92528704735878, 0.0008872371014406283, 0.011535230146151562, -0.012624010212993454, -0.11565111131793405, -0.013764107184334032, -1.031582004225383, -0.0062967862401016655], gphase = [-0.8110030018693761, -3.1098275309525545, 0.837506722408347, -0.8737333121435851, -3.054504399239892, 0.6154931064787439, -0.9234050249748407, -2.9795279373276555, 0.42399923405798245, -0.9745564412063065  …  -0.5827637751656793, -0.4280803614553966, 0.7945337200478205, -0.6584076110384696, 2.571630208733882, -0.5270513733847624, -0.5880975682071283, 0.8214616189586754, -0.8834424963147622, 2.596303149990953]), (c = (params = [1.0430153599579801 1.0838720429706725 … 0.689462552571418 1.5244344518763442; -0.2877845271116076 0.3545628639630058 … 1.090692557567136 1.8751170602620075; … ; 1.62026661473852 1.0787297837415528 … 0.8534911265046562 0.4932316119044106; 0.8396147897316825 1.126965126756855 … 0.8083321402763357 0.7558112062529202], hyperparams = 61.30005852510286), fg = 0.3045751578025416, σimg = 0.9551536143339875, lgamp = [-0.01533351455625373, 0.009748599818479895, -0.9347079523237788, -0.04017504862059126, -0.00037264382543703494, 0.000592470601434628, -0.41390846409251425, 0.1051008877416068, -0.029734803119969944, 0.026441751519692525  …  -0.07165189530607093, 0.006909402635269201, -0.9540201476321472, -0.008709652408477396, 0.005169225861667522, -0.007450121709490281, -0.11668626659335564, -0.0032405598740336146, -1.028578445912532, 0.003504589728335466], gphase = [-0.8090601833786983, 3.096536183735953, 0.8578635407147085, -0.8767542771513969, -3.070187097490882, 0.6578879482197604, -0.933849209882661, -2.9985706677880875, 0.44438618844227407, -0.9807006075444616  …  -0.5833613821907617, -0.4126815511545127, 0.8652307919463053, -0.6464447780734381, 2.628477338836042, -0.5197774714154969, -0.5771724326074652, 0.9070383376048313, -0.8817706609680868, 2.7009899946456035]), (c = (params = [0.6298984401995951 1.4585573187391714 … 1.2042522324752227 0.8945791365874917; 1.527067179357318 0.8677516564810586 … 1.9809136842439767 1.0356932726462922; … ; 1.1435596439647056 1.3698875133911497 … 0.9739030978727546 1.2538755415214953; 1.0791339958298656 0.8965139908038594 … 0.5315184216378409 0.43139567516724664], hyperparams = 55.23222456881166), fg = 0.32220699478692466, σimg = 1.009097094199292, lgamp = [-0.0067747393486831085, 0.013582660151897226, -0.9407093858696816, -0.032322380324249174, -0.0063543038282777775, 0.00274553290850039, -0.33309736948636337, 0.10460023599390565, -0.021412326648546486, 0.019155838606581745  …  -0.08059029333978926, -0.00527574383951005, -0.9386445996019385, -0.004059496925981274, 0.01683706771455022, -0.01701101045633338, -0.1276538148664993, -0.013572168398158595, -1.0332786320515044, -0.010956306655615639], gphase = [-0.8180327360246812, -3.1187144944583505, 0.9231082446672327, -0.8670816209207909, -3.011150081134314, 0.7173976416144422, -0.9300855855631672, -2.941599347577571, 0.5163393052966381, -0.9779428126098789  …  -0.5844458215415512, -0.38086651022219414, 0.8039457049584693, -0.6198640179529274, 2.577566817673359, -0.521432636083851, -0.5288688707081136, 0.810932479297886, -0.8516413855061025, 2.601123874822985]), (c = (params = [-0.03702955238630584 -0.1019484136216008 … 0.2660634973417016 0.7770996660064496; -0.019490566337668554 -0.12563665948229008 … -1.0763656158584247 -0.4798016145617032; … ; 1.117407330986981 0.0479481473399009 … 0.6203100302298363 0.8804104731408329; 0.5009086851005009 0.9968637430478752 … 1.5083851209889323 0.6862837030522836], hyperparams = 12.75253686937719), fg = 0.3087944652346079, σimg = 0.9708916619637458, lgamp = [-0.021909964743065096, 0.02951938633986347, -0.9394917806488077, 0.013799763503329452, 0.009772968378700175, -0.008868628333155543, -0.35970422479419806, 0.11748386242272713, -0.014882332168727678, 0.015293813051770567  …  -0.09212909220953736, 0.0011366298718709965, -0.9531333401765414, -0.006743520321369363, -0.03227719752936948, 0.029704339875466064, -0.07607817436403938, 0.014543587900987823, -1.0039108628934874, 0.01079385615415651], gphase = [-0.8184657856107577, -3.050505470318185, 0.9360855557083304, -0.8672992091964482, -2.9964377065879773, 0.7483131665104344, -0.9246111235584642, -2.913379804317162, 0.5659148330231852, -0.9739585829759938  …  -0.5853981496789975, -0.37516424989592734, 0.7252665410106095, -0.5955146637287493, 2.501887603281616, -0.5211119130892665, -0.5366045814064467, 0.7508306338444077, -0.8092406188750215, 2.538588719794414]), (c = (params = [0.018036876632701194 -0.052267728366586406 … 0.7700360785426872 -0.019116403898442005; -0.5039588167092927 -1.0170737320680099 … 1.1438370332594425 0.15690758012344969; … ; 1.8944626806224925 1.251147130181852 … 1.507416767664301 2.257695946131935; 0.9271595139606933 -0.05078249409062566 … -0.0773447208477155 1.1089882104838842], hyperparams = 19.410047158732922), fg = 0.3149981506463674, σimg = 0.9532560537861618, lgamp = [-0.062751231322591, 0.05590014048515919, -0.9629249455279955, 0.021398623437082417, -0.011034125118740072, 0.01045793722896658, -0.3932824551974471, 0.11689330065374139, 0.016120465910018658, -0.015801057493667003  …  -0.09841982511678181, -0.004095663317112257, -0.958128949906042, -0.023095094116900923, -0.0017309557261940374, -0.0017791814283692115, -0.11036398189843426, -0.004154932017762765, -1.021472210315238, -0.020569140837835592], gphase = [-0.8137085332699231, -3.0591235458010964, 0.9102635274354312, -0.8713981729202677, -2.949302209664687, 0.713736077277578, -0.9260422077030321, -2.877454783721883, 0.5163893341977313, -0.9779983467465686  …  -0.5857088734159023, -0.3543378172857108, 0.7192302318673619, -0.5810026273482225, 2.5025582748043713, -0.5301094884031815, -0.5004206415160435, 0.7712244601676146, -0.7933867396556222, 2.5486017247000934]), (c = (params = [0.27594152740736916 -0.5787321437707614 … 0.8665175529595096 0.21112350834288576; -0.6382964782594133 -1.326639271295839 … 0.9265662549877912 0.1615662953252691; … ; 2.030732417221604 0.9986377982019933 … 1.533191515096395 2.4020061953836214; 0.8885138932981308 0.4192387055373466 … 0.11301260747046724 1.0944988169971073], hyperparams = 17.532732399746926), fg = 0.310716093978145, σimg = 0.9624872621408307, lgamp = [-0.05316614194071191, 0.06038463008875371, -0.9528446782492098, 0.03256303678815633, -0.015506193714024506, 0.018078292683385365, -0.4044987998243128, 0.12797653716691396, 0.009671066516803507, -0.01066522208643972  …  -0.09194174295942405, -0.0049302241854731945, -0.9311485444611565, -0.025821238526801512, -0.0014421652582178124, 0.0020065447528201134, -0.10603355314681052, -0.0019087455414863188, -1.0140882635314656, -0.01587886741451013], gphase = [-0.8128458519058746, -3.0384457463761287, 0.9137187888998938, -0.8716032303927954, -2.979242988619189, 0.7367813069490414, -0.9269102776522129, -2.8833859529427865, 0.5452835191385783, -0.9757804260005412  …  -0.5887708132875814, -0.37163034311652055, 0.7475119575444691, -0.5930202927063833, 2.5248649911983887, -0.5266110497483149, -0.5009434030979162, 0.7944091088241162, -0.8302914268609723, 2.5750833814818943]), (c = (params = [-0.27212037874931977 -0.6678052440807544 … 0.611058183853557 0.9411881317021848; 0.36998162144042734 0.392846735434795 … 0.18003576985456476 0.08138947402123081; … ; 0.9591466589023352 0.42789631869604755 … 0.34942015270279536 1.5121100591301855; 1.4042673742764344 0.039441602960778754 … 0.8323715439873112 1.435371140373917], hyperparams = 28.036361501265276), fg = 0.29620635362835485, σimg = 0.9633979043198713, lgamp = [-0.055082185747768755, 0.054688111663480296, -0.9316845476707899, 0.02904169475700907, -0.014318745365251817, 0.014341513474993678, -0.41442637520585013, 0.1367967823715441, 0.014450800072621782, -0.00965030203723813  …  -0.07709806773154022, 0.004909263408175223, -0.961915110349135, -0.010747020926648148, 0.008492318895986747, -0.012160556842416636, -0.12289637165093671, -0.01134479073380468, -1.0121575461118097, -0.01103250498614686], gphase = [-0.8136538497789378, -2.9875513621909335, 0.906959998098852, -0.8735146006912313, -2.9238461354173815, 0.7073915806006679, -0.9300063938371232, -2.852785789679892, 0.5275706340335099, -0.9755374477116447  …  -0.5846507574596574, -0.3393850618335718, 0.8014479162068159, -0.5394686999264549, 2.570280378430395, -0.5274657604448552, -0.475924615523362, 0.7951376645585205, -0.7682570983396505, 2.5903035540973014])], NamedTuple{(:n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :max_hamiltonian_energy_error, :tree_depth, :numerical_error, :step_size, :nom_step_size, :is_adapt), Tuple{Int64, Bool, Float64, Float64, Float64, Float64, Float64, Int64, Bool, Float64, Float64, Bool}}[(n_steps = 1023, is_accept = 1, acceptance_rate = 0.9926559875781933, log_density = 2936.8356531279674, hamiltonian_energy = -2385.754441047356, hamiltonian_energy_error = 0.008003805291536992, max_hamiltonian_energy_error = 0.01082524344428748, tree_depth = 10, numerical_error = 0, step_size = 0.0001, nom_step_size = 0.0001, is_adapt = 1), (n_steps = 2, is_accept = 1, acceptance_rate = 7.0080726234797135e-37, log_density = 2936.8356531279674, hamiltonian_energy = -2282.903993336208, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 5085.991253484927, tree_depth = 1, numerical_error = 1, step_size = 0.0015545602445261016, nom_step_size = 0.0015545602445261016, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9444658867877068, log_density = 2518.5558121078257, hamiltonian_energy = -2242.1772157051937, hamiltonian_energy_error = 0.06907471067961524, max_hamiltonian_energy_error = 0.11628441818038482, tree_depth = 10, numerical_error = 0, step_size = 0.0003024554438978711, nom_step_size = 0.0003024554438978711, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9385780931336954, log_density = 2271.1136586303523, hamiltonian_energy = -1844.6469250178168, hamiltonian_energy_error = 0.07670650461636797, max_hamiltonian_energy_error = 0.18896597959110295, tree_depth = 10, numerical_error = 0, step_size = 0.0004344283181627208, nom_step_size = 0.0004344283181627208, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9450172483876977, log_density = 2044.9289037388132, hamiltonian_energy = -1636.7318942051375, hamiltonian_energy_error = 0.01278257715057407, max_hamiltonian_energy_error = 0.31878118874647043, tree_depth = 10, numerical_error = 0, step_size = 0.0007010732176278496, nom_step_size = 0.0007010732176278496, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.567278263613081, log_density = 1936.0004689185112, hamiltonian_energy = -1388.273204234487, hamiltonian_energy_error = 0.8645506615009708, max_hamiltonian_energy_error = 7.621851336096142, tree_depth = 10, numerical_error = 0, step_size = 0.0012347123185972456, nom_step_size = 0.0012347123185972456, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.756791993387574, log_density = 1742.9576897668633, hamiltonian_energy = -1209.023130861764, hamiltonian_energy_error = -0.14547307695283962, max_hamiltonian_energy_error = 1.0538280520240733, tree_depth = 10, numerical_error = 0, step_size = 0.0007096790842531722, nom_step_size = 0.0007096790842531722, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9828447454733104, log_density = 1557.2578905691394, hamiltonian_energy = -1066.8867968380682, hamiltonian_energy_error = -0.02895338243229162, max_hamiltonian_energy_error = -0.42719310475627026, tree_depth = 10, numerical_error = 0, step_size = 0.0007207307685665362, nom_step_size = 0.0007207307685665362, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.7317486358065133, log_density = 1467.9668785820727, hamiltonian_energy = -916.332269562916, hamiltonian_energy_error = 0.4059567187044877, max_hamiltonian_energy_error = 1.872758396074346, tree_depth = 10, numerical_error = 0, step_size = 0.0014934805939154647, nom_step_size = 0.0014934805939154647, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.41955351645518524, log_density = 1330.0979440069432, hamiltonian_energy = -788.4586406601729, hamiltonian_energy_error = -1.4350000730820511, max_hamiltonian_energy_error = 12.391593013903162, tree_depth = 10, numerical_error = 0, step_size = 0.0014125735079172012, nom_step_size = 0.0014125735079172012, is_adapt = 1)  …  (n_steps = 511, is_accept = 1, acceptance_rate = 0.913475040373775, log_density = 1141.6183199218374, hamiltonian_energy = -471.39973551799403, hamiltonian_energy_error = 0.24549929666363823, max_hamiltonian_energy_error = 0.6713993175856103, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.8585208930107143, log_density = 1126.2793407720776, hamiltonian_energy = -428.94735706509186, hamiltonian_energy_error = -0.07435567476727556, max_hamiltonian_energy_error = 0.7982720563570638, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.972952145813327, log_density = 1137.8104365715042, hamiltonian_energy = -433.35031739960084, hamiltonian_energy_error = -0.13042783656158008, max_hamiltonian_energy_error = -0.5834098566583634, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.8360557060942708, log_density = 1140.644714978218, hamiltonian_energy = -442.8918333005423, hamiltonian_energy_error = 8.001539879387565e-5, max_hamiltonian_energy_error = 0.677833574259239, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9524293260815583, log_density = 1142.828554757725, hamiltonian_energy = -475.714053798599, hamiltonian_energy_error = 0.09052517893246659, max_hamiltonian_energy_error = -0.5181884322630594, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.6397541890860913, log_density = 1188.3858683513665, hamiltonian_energy = -473.0749038349974, hamiltonian_energy_error = -0.01644521444734437, max_hamiltonian_energy_error = 1.1279425268362502, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9574734962924285, log_density = 1173.2622134297496, hamiltonian_energy = -490.9700366987913, hamiltonian_energy_error = 0.10124278350656368, max_hamiltonian_energy_error = -0.4390620710021267, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.6474247264562024, log_density = 1165.5337010289952, hamiltonian_energy = -525.1016638435669, hamiltonian_energy_error = 0.6824202513382716, max_hamiltonian_energy_error = 1.4039411507762907, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9403534043381111, log_density = 1162.287232896325, hamiltonian_energy = -513.6326674014422, hamiltonian_energy_error = 0.12546669554160417, max_hamiltonian_energy_error = -1.0832590268414606, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.8651348580406518, log_density = 1143.1534289770245, hamiltonian_energy = -474.01300804794687, hamiltonian_energy_error = -0.8171702271112053, max_hamiltonian_energy_error = -1.279603443184783, tree_depth = 9, numerical_error = 0, step_size = 0.005854827776785766, nom_step_size = 0.005854827776785766, is_adapt = 0)])
Note

The above sampler will store the samples in memory, i.e. RAM. For large models this can lead to out-of-memory issues. To fix that you can include the keyword argument saveto = DiskStore() which periodically saves the samples to disk limiting memory useage. You can load the chain using load_table(diskout) where diskout is the object returned from sample. For more information please see ComradeAHMC.

Now we prune the adaptation phase

chain = chain[501:end]
Table with 5 columns and 200 rows:
      c                     fg        σimg      lgamp                 ⋯
    ┌──────────────────────────────────────────────────────────────────
 1  │ (params = [-0.36326…  0.200902  0.860534  [-0.00356013, -0.00…  ⋯
 2  │ (params = [0.139049…  0.211479  0.875864  [-0.00803269, 0.011…  ⋯
 3  │ (params = [0.820433…  0.199472  0.857334  [0.00564365, -0.003…  ⋯
 4  │ (params = [0.555763…  0.208868  0.875653  [0.00186956, 0.0012…  ⋯
 5  │ (params = [0.546449…  0.20403   0.855899  [0.000134808, -0.00…  ⋯
 6  │ (params = [0.044505…  0.219803  0.930671  [-0.00638664, 0.003…  ⋯
 7  │ (params = [0.270818…  0.204429  0.908768  [-0.0211515, 0.0193…  ⋯
 8  │ (params = [0.171623…  0.212662  0.912665  [-0.0169229, 0.0185…  ⋯
 9  │ (params = [1.28162 …  0.213062  0.95642   [-0.00543214, 0.008…  ⋯
 10 │ (params = [0.921775…  0.21592   0.908115  [-0.0135997, 0.0160…  ⋯
 11 │ (params = [0.776145…  0.223757  0.944354  [-0.00624355, 0.005…  ⋯
 12 │ (params = [0.569273…  0.211688  0.945337  [0.00112532, 0.0003…  ⋯
 13 │ (params = [1.16881 …  0.221181  0.90916   [-3.9746e-5, 0.0011…  ⋯
 14 │ (params = [2.21913 …  0.22024   0.990129  [0.0132194, -0.0103…  ⋯
 15 │ (params = [0.892848…  0.245825  1.0278    [0.00171694, 0.0067…  ⋯
 16 │ (params = [1.45334 …  0.264395  1.02395   [-0.0154122, 0.0116…  ⋯
 17 │ (params = [1.45822 …  0.262268  0.98737   [-0.0207525, 0.0221…  ⋯
 ⋮  │          ⋮               ⋮         ⋮               ⋮            ⋱
Warning

This should be run for likely an order of magnitude more steps to properly estimate expectations of the posterior

Now that we have our posterior, we can put error bars on all of our plots above. Let's start by finding the mean and standard deviation of the gain phases

gphase  = hcat(chain.gphase...)
mgphase = mean(gphase, dims=2)
sgphase = std(gphase, dims=2)
104×1 Matrix{Float64}:
 0.0036607480404510346
 0.9572636458060648
 0.2387292007473356
 0.0030414295979255394
 0.1713598278136156
 0.24133820845764714
 0.0031935588408206843
 0.17550783401033032
 0.24253086207062938
 0.003229558214215475
 ⋮
 0.22287383536955516
 0.30576653586926567
 0.17602426517044023
 1.6777253271501564
 0.003789397604349588
 0.2183800319789883
 0.30134028816625275
 0.17187015516602824
 1.7962465669942818

and now the gain amplitudes

gamp  = exp.(hcat(chain.lgamp...))
mgamp = mean(gamp, dims=2)
sgamp = std(gamp, dims=2)
129×1 Matrix{Float64}:
 0.018831672803544076
 0.01919568886631146
 0.030591740354320957
 0.041444692448386555
 0.01592645700996742
 0.01619896017686868
 0.0390527700141429
 0.04342275073075323
 0.0160808024699145
 0.016733885101951146
 ⋮
 0.008737772012239527
 0.017942399858837294
 0.009532788796534325
 0.014858671260356065
 0.0156580270931832
 0.023286209696402906
 0.009504994320670535
 0.018137467162359536
 0.009184976167852394

Now we can use the measurements package to automatically plot everything with error bars. First we create a caltable the same way but making sure all of our variables have errors attached to them.

using Measurements
gmeas_am = measurement.(mgamp, sgamp)
ctable_am = caltable(gcache, vec(gmeas_am)) # caltable expects gmeas_am to be a Vector
gmeas_ph = measurement.(mgphase, sgphase)
ctable_ph = caltable(gcachep, vec(gmeas_ph))
───────────┬────────────────────────────────────────────────────────────────────
      time │      AA            AP          AZ           JC           LM       ⋯
───────────┼────────────────────────────────────────────────────────────────────
 0.917±0.0 │ 1.0±0.0  -0.814±0.004     missing      missing   -2.75±0.96   0.9 ⋯
 1.217±0.0 │ 1.0±0.0  -0.871±0.003     missing      missing   -2.82±0.17   0.7 ⋯
 1.517±0.0 │ 1.0±0.0  -0.929±0.003     missing      missing   -2.75±0.18   0.5 ⋯
 1.817±0.0 │ 1.0±0.0  -0.977±0.003     missing      missing    -2.7±0.18   0.3 ⋯
 2.117±0.0 │ 1.0±0.0  -1.002±0.005     missing      missing   -2.63±0.18   0.1 ⋯
  2.45±0.0 │ missing       1.0±0.0     missing      missing      1.6±2.4    1. ⋯
  2.75±0.0 │ 1.0±0.0  -1.045±0.007     missing      missing   -2.56±0.19   -0. ⋯
  3.05±0.0 │ 1.0±0.0  -1.132±0.003     missing      missing   -2.53±0.19  -0.3 ⋯
  3.35±0.0 │ 1.0±0.0  -1.152±0.003     missing      missing    -2.54±0.2  -0.5 ⋯
 3.683±0.0 │ 1.0±0.0  -1.164±0.005     1.8±2.2      missing   -2.59±0.19  -0.7 ⋯
 3.983±0.0 │ 1.0±0.0  -1.167±0.003     2.2±1.7      missing    -2.64±0.2  -0.9 ⋯
 4.283±0.0 │ 1.0±0.0   -1.16±0.004   2.59±0.23      missing    -2.71±0.2  -1.0 ⋯
 4.583±0.0 │ 1.0±0.0  -1.143±0.003   2.41±0.23    -0.51±0.3   -2.75±0.62  -1.2 ⋯
 4.917±0.0 │ 1.0±0.0  -1.073±0.007    2.2±0.23    -0.26±0.3     -1.8±2.3  -1.4 ⋯
 5.183±0.0 │ 1.0±0.0  -1.078±0.005   2.04±0.23  -0.092±0.31     0.16±3.0  -1.6 ⋯
  5.45±0.0 │ 1.0±0.0  -1.116±0.012    1.9±0.24    0.08±0.32      1.9±2.2  -1.7 ⋯
     ⋮     │    ⋮          ⋮            ⋮            ⋮            ⋮            ⋱
───────────┴────────────────────────────────────────────────────────────────────
                                                    2 columns and 9 rows omitted

Now let's plot the phase curves

plot(ctable_ph, layout=(3,3), size=(600,500))
Example block output

and now the amplitude curves

plot(ctable_am, layout=(3,3), size=(600,500))
Example block output

Finally let's construct some representative image reconstructions.

samples = skymodel.(Ref(post), chain[begin:2:end])
imgs = intensitymap.(samples, fovx, fovy, 128,  128)

mimg = mean(imgs)
simg = std(imgs)
fig = CM.Figure(;resolution=(400, 400))
CM.image(fig[1,1], mimg,
                   axis=(xreversed=true, aspect=1, title="Mean Image"),
                   colormap=:afmhot)
CM.image(fig[1,2], simg./(max.(mimg, 1e-5)),
                   axis=(xreversed=true, aspect=1, title="1/SNR",),
                   colormap=:afmhot)
CM.image(fig[2,1], imgs[1],
                   axis=(xreversed=true, aspect=1,title="Draw 1"),
                   colormap=:afmhot)
CM.image(fig[2,2], imgs[end],
                   axis=(xreversed=true, aspect=1,title="Draw 2"),
                   colormap=:afmhot)
fig

Now let's check the residuals

p = plot();
for s in sample(chain, 10)
    residual!(p, vlbimodel(post, s), dvis)
end
p
Example block output

And viola, you have just finished making a preliminary image and instrument model reconstruction. In reality, you should run the sample step for many more MCMC steps to get a reliable estimate for the reconstructed image and instrument model parameters.

This page was generated using Literate.jl.